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Enthalpy Damping for the Steady Euler Equations 

DENNIS C. JESPERSEN 1 

MS 202A-1, NASA Ames Research Center 
Moffett Field, CA 94035 USA 

SUMMARY. For inviscid steady flow problems where the enthalpy is constant at steady state, it has been proposed by Jameson, 
Schmidt, and Turkelto use the difference between the local enthalpy and the steady state enthalpy as a driving term to accelerate 
convergence of iterative schemes This idea is analyzed here, both on the level of the partial differential equation and on the 
level of a particular finite difference scheme It is shown that for the two-dimensional unsteady Euler equations, a hyperbolic 
system with eigenvalues on the imaginary axis, there is no enthalpy damping strategy which can move all the eigenvalues into 
the open left half plane. For the numerical scheme, however, the analysis shows and examples verify that enthalpy damping can 
be effective in accelerating convergence to steady state. 


l. Introduction 

Rapidly convergent numerical schemes for steady state solutions of the inviscid compressible fluid dynamics equa- 
tions, the Euler equations, have been studied mtensively in the last few years. Algorithms and codes have been 
developed by Jameson, Schmidt, and Turkel (non-multigrid) (Ref 1); Jameson (multigrid) (Ref 2), Rizzi (Ref 3); 
Pulliam (Ref. 4), Lerat, Sides, and Daru (Ref 5), Johnson (Ref. 6); Ni (Ref. 7); Agarwal and Deese (Ref 8); Osher 
and Chakravarthy (Ref. 9), Mulder and van Leer (Ref 10), and by Bunmg and Steger (Ref 11), among others One 
of the most successful and widely used algorithms is embodied in the FL052 codes of A Jameson The algorithm 
used for the Euler equations in these codes is based on an explicit multistage time-stepping scheme, with (for steady 
state problems) acceleration devices such as local time step selection, implicit residual averaging, multignd, and en- 
thalpy damping It is the purpose of this paper to provide an analysis of enthalpy damping for the two-dimensional 
Euler equations 

The original derivation of enthalpy damping by Jameson, Schmidt, and Turkel (Ref. 1) was based on reducing the 
unsteady Euler equations to the unsteady potential equation (a wave equation), introducing frictional damping to 
modify the equation to the telegraph equation, and translating the new equation back to the Euler equations The 
next section begins by reviewing this derivation It can be verified that enthalpy damping is effective in enhancing 
convergence for the Euler equations, and an example calculation is included. 

The two-dimensional Euler equations with enthalpy damping are studied in section 3. The Euler equations are a 
system of partial differential equations of the form 

w t + dj£(w) -I- d y J(w) -f a)i(w) =0 (11) 

where the enthalpy damping is embodied in the undifferentiated term W, and a determines the strength of the 
damping term We linearize (1 1) around a state tv 0 with constant enthalpy, obtaining a linear variable coefficient 
system 

wt + d x (A(w 0 (z, y))w) + d v (B(w 0 (x, y))w) + aC(w 0 (z, y))w = 0 (12) 

where tv is the perturbation from w 0 and A = d£/dtv , B — d? fdw, C = dMfdtv . We study the frozen coefficient 
problem 

tv t + Adjtv + Bdytv + aCw = 0 (13) 

Fourier transforming (1 3) gives the system 

^ + t(uiA + u 2 B)w + aCw = 0 (14) 

where w(x, y, f) = e *( w i*+ w 3y)u)(f]. For a = 0, the system (1 3) is hyperbolic and all solutions tv of (1 4) are of the 
form w = Ylj e Aj< «^(0) where A is purely imaginary. For a > 0 the natural conjecture is that the enthalpy damping 
helps by pulling the A ; into the left half of the complex plane. We show that of the four eigenvalues A_, for any 
pair (^ 1 ,^ 2 ), one eigenvalue always remains on the imaginary axis, so the system (1 3) does not have the property 
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that all its solutions decay to zero. The analysis also suggests a form of the enthalpy damping term for the energy 
equation more natural than that of Jameson, Schmidt, and Turkel (Ref. 1). 

If the Euler equations at the partial differential equation level do not have decaying solutions when enthalpy 
damping is in effect, why do the numerical results show that enthalpy damping is useful? This question is addressed 
in section 4. To get a flavor of the idea, consider the system of ordinary differential equations 

du , 

— + Au = 0 (15) 

where the complex matrix A is assumed to be diagonalizable with all its eigenvalues on the imaginary axis, say at 
1 A 1 , 1 A 2 , .. , with 0 < Ax < A 2 < ... . Solutions to (1 5) do not decay to zero as t — ► 00 . However, if (15) is 
integrated with the implicit Euler method 


«n+i - « ft + A/A« ft+l = 0 (1.6) 

the numerical solution ti n will go to zero as n — ♦ 00 because the eigenvalues of A are contained in the interior of the 
stability region of the method Furthermore, the eigenvector associated with the eigenvalue i\ x will decay slowest 
because its damping factor is closest to 1 Thus if the matrix -4 is somehow modified so that the eigenvalue Aj moves 
to the left, such that the modified eigenvalue is farther from 1 in the complex plane than the original eigenvalue, and 
such that all the other eigenvalues are unchanged, then the numerical solution of the modified system will decay to 
zero faster than the numerical solution of the original system In section 4 we attempt to show that roughly the same 
phenomenon occurs with the Euler equations and an explicit multistage integration scheme (1 e , the eigenvectors 
obstructing the convergence are associated with eigenvalues which are moved to the left when enthalpy damping is 
used) Thus, even though one of the eigenvalues does not move off the imaginary axis, the convergence rate can 
still improve An attempt is made to quantify the improvement in convergence rate and to find an optimal damping 
parameter a 

The major omission of this work is the neglect of numerical dissipation. The numerical scheme is actually an 
approximation not to (1 1), but to 


tv t + d x £(w) + d y 7(w) + aJ/(u/) + cD(w) = 0 (1.7) 

where D is a dissipation operator (e g , a blend of second- and fourth-order dissipation terms) The inclusion of a 
dissipation term makes the analysis much more cumbersome, a point which will be touched on again If an analysis 
including the dissipation term could be carried out, results more relevant to the performance of the computer codes 
could be obtained 

The author has had the benefit of several illuminating discussions with Eli Turkel in the course of this work Many 
of the algebraic computations in sections 3 and 4 were carried out with vaztma, the VAX version of Macsyma (Ref. 
12) Some of the computations would have been extremely difficult or impossible without vazxma. 


2. Derivation of Enthalpy Damping 

This section reviews the derivation of enthalpy damping via the unsteady potential equation and telegraph equation 
as originally sketched in Jameson, Schmidt, and Turkel (Ref 1) The two-dimensional Euler equations in conservation 
form for a perfect gas may be written 


Pt + (p«)i + (pv) v = 0 
{pu) t + (pu 2 + p) x + (puv) y =0 
(pv)t + (put;)* + (pt; 2 +p)„ = 0 
(pE) t + (puH)j + (pvH)y = 0 

where p is the density, u and v are Cartesian velocity components, E is total energy per unit mass, p is pressure, 
and H is the enthalpy, defined by 

H = E + pfp 
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One can verify that 


D_ n= 1£P 
Dt p dt 

where D/Dt is the material derivative D/Dt := d/dt + v • V. (The notation a •= b means a is defined as b ) Hence 
in steady flow H is constant along streamlines, and if flow conditions are uniform at upstream infinity, then H is 
everywhere constant at steady state. 

In nonconservative form the Euler equations may be written 


Pi + up* + vp v + p(« i + v y ) = 0 
«{ + uu x + vu y + p x /p = 0 

vt + uv x + vv y + p y /p = 0 (2 2) 

„ 1 

Et + UE X + VEy + -(« X + Vy) + “(«p x + VPy) SS 0 

P P 

The fourth equation here can be replaced by the transport equation for entropy 

St 4* uS x 4* vSy = 0 

The choice of the fourth equation affects the form of enthalpy damping that is derived 

If the flow is barotropic (i e., p = p(p)) one can show by differentiating the i-momentum equation with respect 
to y, the y-momentum equation with respect to x, and subtracting the latter from the former, that the vorticity 
u = v x — u y satisfies 

cjf + (wu) x 4* (wi?)y = 0 

With this and the use of the continuity equation one can show that 

+ ud 2 (u/p) + vd v (ujp) = 0 (2 3) 

Thus if the fluid is irrotational (i e , u; = 0) at time t — 0 then it is irrotational for all t > 0. For an irrotational 
fluid a velocity potential (j> may be introduced, where <j> x = u and <p y = v Then the momentum equations may be 
integrated to give the unsteady Bernoulli equation 


dt 




+ 


/?= G « 


(2 4) 


where G is an arbitrary function of t If the flow is homentropic (entropy everywhere constant), then we have p = Ap 1 
and / dp/p = 7 p /[(7 — l)p] For an ideal gas we have p = (7 - 1 )p\E — (n 2 4 - v 2 )/2] where 7 is the ratio of specific 
heats Hence 


r£ + 


u 2 4 * V 2 
2 


and the unsteady Bernoulli equation becomes 


d<}> 

dt 


+ H = G(/) 


For exterior flow problems, if the conditions are uniform at infinity, this becomes 




(2 5) 


where H is the enthalpy at infinity 
An alternate form of the unsteady Bernoulli equation (2 5) is 


. + 

4>t + — ^ — 


7-1 


Ap - 1 - 1 = H 0 
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To obtain the unsteady potential equation, this equation is differentiated with respect to t , x, and y, and solved 
for pt , p 2 , p y . These expressions for p*, p*, and p y are substituted into the mass conservation equation. Using the 
relation giving the speed of sound as c 2 = 7 p/p, one obtains the unsteady potential equation 

(f>tt + 2 u<t> xi + 2r0 y « = (c 2 - u 2 )^* - 2uv^ v + (c 2 - t? 2 )0 yy 

For constant u and t; the change of coordinates £ = x — tit, 1/ = y — t?f transforms this to the wave equation 

<l>tt = c 2 (^e 

A related equation (Ref. 13) is the telegraph equation 

<l>tt + a<t>t = c 2 (<^ 4- (f>Tjn) 

For a > 0 the solutions of this equation decay as t 00. Transforming this equation back via x = £ + uf , y = rj 4- t?J 
we obtain the damped version of (2 6), 

4>tt 4- 2 u<j> xt 4- 2 v(j> yt 4* a<£* = (c 2 - u 2 )^ - 2uv<t> xy 4- (c 2 - t? 2 )0 yy (2.7) 

This is the equation one would have obtained from the mass conservation equation if the mass conservation equation 

had been 

pt 4- up* 4- vpy 4 - p{u x 4 - Vy) = ap<pt 

Now, (j> is unknown, but the unsteady Bernoulli equation (2 5) allows us to replace (f> t by Hoc - H, obtaining the 
modified mass conservation equation 

Pt 4- up x 4- vp v 4- p(u x 4- r y ) 4- ap(H - H 00 ) = 0 (28) 

Using this in place of the first equation of (2 2) and recombining into conservation form, one obtains the system of 

equations 

Pt 4- (pu) x 4- (pu)y 4- otp(H - 
(pu)f 4- (pu 2 4* p)j 4- (puv) y 4- apu(H - 
(pt?)* 4- (puv) x 4* (pt? 2 +p) y 4- apt?(if - 
(pE) t 4- (puH) x 4* (pvH)y 4- apE(H - 

If the entropy equation is used as the fourth equation of (2 2), recombining 
fourth equation in (2 9) Then (2 9) is replaced by 

Pt + (pu)* 4- (pv) y 4- ap(H - 
(P«h 4- (pu 2 4* p) x 4- (p«t?) y 4- apu(H - 
(pt?)* 4- (puv) x 4- (pt? 2 4- p)y 4- apv(H - 
(pE) t 4- (puH) x 4- (pvH)y 4- apH(H - 

In the next section we will analyze a general form of enthalpy damping which includes both (2 9) and (2 10). 

In summary, for subsonic irrotational flow, systems (2 9) and (2 10) give rise (locally) to the telegraph equation 
while the system (2 1) gives rise (locally) to the wave equation. Furthermore, we may hope that steady states of 
(2 9) and (2 10) are steady states of (2 1), and vice versa. It is clear that any steady state of (2 1) is also a steady 
state of (2.9) and (2 10) (since steady states of (2 1) have constant enthalpy). It is not clear that steady states of 
(2 9) are steady states of (2 1). It is true, however, that steady solutions of (2 10) have enthalpy constant along 
streamlines, hence if flow conditions are uniform at upstream infinity then steady solutions of (2 10) have enthalpy 
constant everywhere, hence steady solutions of (2 10) are also steady solutions of (2 1). To see that steady solutions 
of (2 10) have enthalpy constant along streamlines, multiply the first equation of (2 10) by H and subtract the result 
from the fourth equation of (2 10). If the time derivatives are zero one finds that uH x 4- vH y = 0, hence H is constant 
along streamlines 


(2 9) 


H 00) =0 

ffoo) =0 
Hoo) =0 
ffoo) =0 

into conservation form produces a different 


ffoo) =0 
ffoo) =0 
ffoo) = 0 
ffoo)=0 


(2 10 ) 


( 2 . 6 ) 
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In place of the fourth equation of (2.9) we can use 

(p(E - ffoo))i + (pti(J5T - JToo)), + (pv(H - ffoo)) y 4- <*p(E - H^H - Hoo) = 0 (2.11) 

which is obtained by multiplying the first equation of (2.9) by H oo and subtracting from the last equation of (2.9). 

It is suggested by Jameson, Schmidt, and Turkel (Ref. 1) that a forcing term apH (H-Hoo) in the fourth equation 
of (2 9) “can be destabilizing,” and it is suggested that this term be replaced with a(J/ - ffoo)- The computer codes 
seen by this author use a forcing term ap(H - /Too) in (2 11); i.e., E - H <*> is replaced by 1 in the fourth term 
of (2 11). This ad hoc fix seems undesirable, if only on dimensional grounds. One of the purposes of this paper is 
to rationally derive a reasonable form for the forcing term in the energy equation, thereby clearing up some of the 
confusion that seems to exist. 

Figure 2 1 shows a comparison of a Euler equation calculation with and without enthalpy damping. The code 
used was FL052R (supplied by A Jameson) with the multigrid and implicit residual averaging options turned off. 
The initial solution was a partially converged solution obtained by gnd continuation (interpolation from a coarser 
grid) The enthalpy damping factor a was 0 05 The convergence rate improved from 0 9966 without enthalpy 
damping to 0 9954 with enthalpy damping. Jameson, Schmidt, and Turkel (Ref. 1) show an example of improvement 
in convergence rate from 0 9909 to 0 9863. The greater gain reported there may perhaps be due to the smaller 
number of grid points (128 x 32) or other mesh differences, differences in the boundary conditions, or differences in 
implementation of the numerical algorithm 


3 Analysis of Enthalpy Damping for the Differential Equation 
Consider the system of partial differential equations 

tv t + d x E (tv) + dy7(w) + aJ/(u?) =0 (3 1) 


where 

w = ( p,pu,pv,pE) T 
£ = ( pu , pu 2 + p, puv, puH) T 
7 = (pv,puv,pv 2 + p,pvH) T 

Jf = ( P(H - Hoo), P«(H - Hoo), pv{H - Hoo), phi(H - 


(3 2) 


Here we have taken a fairly general form of the fourth component of )l , with h i — hi(tv) On dimensional grounds, 
h 4 should have the same units as E 

The first step is to linearize about a steady solution tv 0y and since we are interested m steady solutions with 
constant enthalpy, we will assume Wq has constant enthalpy The first-order terms in the equation for the variation 
t^o + w then give the linear variable coefficient equation for tv 


d t w + d x (A(wo(x , y))w) + d y (£(u> 0 (z, y))w) + aC(w 0 (x,y))w = 0 (33) 

where A, JB, G are the Jacobian matrices of £ , /, and J/, respectively We will study this variable coefficient system 
by freezing the coefficients at an arbitrary (xo,Uq) This gives the constant coefficient problem 

dttv + Ad x w 4- Bdytv 4- aCw — 0, (3 4) 

The justification for this reduction is the conjecture that if all the constant coefficient problems (3 4) ha\e w — ► 0 
as t — * oo, then the variable coefficient problem (3 3) has tv — ► 0 as t — ► oo, and that if the variable coefficient 
problem (3 3) has tv — ► 0 as t — ► oo then the nonlinear problem is (locally) asymptotically stable at a steady state 
with constant enthalpy This author is not aware of any theorems which rigorously quantify these conjectures 
We will analyze (3 4) by Fourier methods, take 

w(x,y,t)=e^ u>x+u ^w(t) 
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Figure 2 1 L 2 norm of residual in density equation for FL052R without 
enthalpy damping (solid line) and with enthalpy damping (dashed line) 
NACA0012 airfoil, A/qo = 0 8, angle of attack= 1 25°, 192 x 32 mesh, 
CFL=2, four-stage method, coefficients (1/4, 1/3, 1/2, 1), two dissipation 
evaluations Iteration started with partially converged solution. 


and substitute in (3 4) to obtain the ordinary differential equation 

^ + [i(u; 1 A + u 2 B) + aC\w = 0 (3 5) 

Here take tv = e zi tD where z is complex and w is a constant four-vector. This gives 


(zl H- i(u/j A + u 2 B) *f aC)tD = 0 (3 6) 

If all the zeros of the polynomial p(z) := det(z/ + i(cj 1 A + u 2 B) + aC) lie in the left half-plane {*( 2 ) < 0}, then all 
solutions of (3 4) decay to 0 as t — ► 00 . 

Note that C has rank one, since we have tf(u>) = (ur t , w 2y tvJi A ) T (H - Hoo) and 


G 


21 

dtv 



U>1 \ 


tv 2 

H=Hoo 

w i 


J 
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The zeros of the polynomial p(z) are -1 times the eigenvalues of the matrix i((jj x A + u 2 B) + aC . We know that 
i{u x A + u 2 B) is diagonalizable with eigenvalues iuU, iuU, iu(U + c), and iu(U - c), where U = u x u + u 2 v and 
u 2 = u\ + The following is a theorem about rank one perturbations of diagonalizable matrices with multiple 
eigenvalues 

THEOREM. Let R be a (complex) matrix with a double eigenvalue A and two linearly independent eigenvectors 
associated with A. Let S have rank one. Then R + S has A as an eigenvalue. 


PROOF Let tq, v 2 be two linearly independent eigenvectors associated with A. Since 5 has rank one, there is some 
nontrivial combination av x 4* bv 2 such that S(av x + bv 2 ) = 0. Then (R + S)(atq + bv 2 ) = A(avi + bv 2 ) 

Applying the theorem with R = t(u x A + u 2 B) and S = aO, we see that t(u x A 4- u 2 B) *f aC has an eigenvalue 
iuU. This leads to the conclusion that enthalpy damping cannot produce a system (3.4) which has solutions that 
decay to 0 as t — ► oo. This result is very general; it applies to any form of enthalpy damping for which the enthalpy 
damping term is proportional to the difference between the local enthalpy and the freestream enthalpy Indeed, if 
the enthalpy damping term in (3 2) is formulated as H = (H - BooM/M 1 ^)* fi 2 (w) y 0s(tv), 0 4 (tv)) T y then the matrix 
G has rank one. 

We are going to further study the system of equations with enthalpy damping, to investigate the behavior of the 
eigenvalues that do move off the imaginary axis and, to find a reasonable choice for h 4 To do this, it is not convenient 
to work with the matrices A and B, for the entries of these matrices are complicated functions of (w\ y tv 2 , tv 4 ) It 
is easier to work with similarity transforms N~ 1 AN and N~ l BN , where N is the matrix dtv/dtv and w = (p, u, v,p). 
This simply amounts to writing the Euler equations in primitive variables rather than in conservative variables One 
has 


i 

u p 

v Op 

^ (u 2 + v 2 )/2 pu pv 


\ 

1/(7 -1)7 


N~ l 


1 

- «/p i Ip 

- v/p 0 

V(7~ l)(u 2 + t; 2 )/2 (l-7)« 


\ 

i Ip 

(l-7)v 7-1 J 


The matrices N~ l AN and N~ l BN have been given m Warming and Beam (Ref. 14) They are 


N~ l AN = 


(» 

p 

0 

0 A 


fv 

0 

p 

0 V 

0 

U 

0 

1 Ip 

, N~ l DN = 

0 

V 

0 

0 

0 

0 

u 

0 

0 

0 

V 

1/p 

VO 

pc 2 

0 

« ) 


U 

0 

pc 2 

v J 


The matrix N~*CN turns out to be 


N~'CN = 


( C V(1 — 7) pu pv 7 /( 7 -!) ^ 

0 0 0 0 

0 0 0 0 

^ < first row times (1 — 7)(u 2 + v 2 - 2h 4 )/2 > J 


Since p{z) = det(zl + i(uj x N~ 1 AN + u 2 N~ l BN) + aN~ l CN), we can use the expressions for N~ X AN, N~ i BN 1 
and N~ l GN to compute p{z). It turns out that with U = u? x u-\-w 2 v and w 2 .= we have p{z) — Po( z )d~oq(z) 

where 

Pot*) = {z + iuU) 2 (z A iu(U + c)){z + iw(f7 - c)) 
q(z) = 7 (hi - E)z s + iwf7((27 + l)(h 4 - E) - c 2 /i)z' 2 

w 2[£l( 2 tf2 _ + (^ 4 ~ — (c a - (7 + 2 )U 2 )\z ( 3 7 ) 

-iw s U(U a -c 2 )(h 4 -E-£) 
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The polynomial po(z) b the polynomial one would obtain with no enthalpy damping; its roots are purely imaginary, 
which shows that the Euler equations without enthalpy damping are a hyperbolic system. (To show hyperbolicity, 
one also has to verify that there are linearly independent eigenvectors associated with the double zero at z = -twU.) 
It can be verified that g(-iu;f7) = 0; this b a consequence of the theorem stated above. 

We can obtain further information. First, take U = cM and replace z by ccjf. Next, take h A - E = kc 2 / 7. Note 
that for k = 0 we have = E and for k = 1 we have ft 4 = H. The polynomial p[z) becomes 

w 4 c 4 [(f + <A/) 2 (f + i(M + l))(f + 1 (A/ - 1)) 

+ — {*f s + — A/(2&7 + k - l)f 2 

(aJ 7 

+ i((2 -2k- *i)A/ 2 +k- l)f - - 1)A/(A/ 2 - 1)}] 


This factors as 

wV(f + iA/){f 8 + 3«Mf 2 + (1 - 3A/ 2 )f + iA/(l - A/ 2 ) + — [fcf 2 + i— (7ft + * - l)f + £— i(l - A/ 2 )]} 

w 7 7 (3 8) 

=: w 4 c 4 (f + iAf)J2(f) 

(The polynomial R($) is exactly the polynomial one would obtam by analyzing enthalpy damping applied to the 
Euler equations m one space dimension.) The plan now b to study the location of the zeros of the polynomial R(f ). 
Assume o > 0 Writing /2(f) = 5(f) + (ac/u)T(s ), we note that for f on the imaginary axis, 5(f) is purely imaginary 
and T(f ) is real, hence R cannot vanish on the imaginary axis unless 5 and T have a common zero on the imaginary 
axis. The zeros of 5 he on the imaginary axb and are -lAf, -i(M + 1), -i(Af - 1) At these f, T(f) takes on the 
values (k - l)/7, -(A/ + 1)(7&- l)/ 7, and (Af - l)(7* “ *+ l)/7> respectively. Thus if A/ 2 ^ 1 and k ^ 1, R has 
no zeros on the imaginary axis Hence in this case we may choose any convenient value of M and use the algorithm 
given by J J. H. Miller (Ref. 15) for deciding whether a polynomial has all its roots in the open left half-plane. 

The algorithm is as follows For any function p(z), define p*(z) = p(-T) Given a polynomial p(z), take a point f 
in the left half-plane for which p[z) ^ 0 ^ p*(z), an( ^ define 

p, z \ _ p*(f)p(*) -p(?)p*( g ) 

* ~ r — f 

Then p(z) is a polynomial of degree less than p[z) The theorem is that p(z) has all its zeros in the open left 
half-plane if and only if |p*(f)l > |p(f)l aQ d P has all its zeros in the open left half-plane Thus one can recursively 
apply this theorem, the degree of the polynomial decreasing at each step 

Consider the case A/ 2 < 1. In thb case M = 0 is a convenient value to pick, and one can show (using Miller’s 
algorithm) that R has all its zeros in (9?(f) < 0} if and only if k > 1. Now consider the case A/ 2 > 1 In this case 
one can choose A/ = 2 (for example), the result is that R never has all its zeros in the open left half-plane Nou 
consider the case A/ 2 = 1 In this case R has a zero at f = 0, and one can show that the other zeros of R are in the 
open left half-plane if and only if k > 1. In summary, the polynomial R($) has all its zeros in the open left half-plane 
(9£(£) < 0} if and only if A/ 2 < 1 and k > 1 Note that the condition A/ 2 < 1 is (cji« + cj 2 u ) 2 < c 2 ( + cj|)t which 
holds for all (wj, oj 2 ) if and only if u 2 + v 2 < c 3 , i e , the flow is subsonic It is interesting to note further that k = 0 
(unstable) is the case h A = E, the formulation of enthalpy damping we obtained in section 2 from the potential 
equation with enthalpy damping, while the case k = 1 (neutrally stable) is the same as h A = H , the formulation that 
“can lead to instability” according to Jameson, Schmidt, and Turkel (Ref 1) 

Thus enthalpy damping leads to a differential system which has a general solution that does not decay to zero 
Of the four eigenvalues on the imaginary axis when there is no enthalpy damping, one stays on the imaginary axis 
when enthalpy damping is introduced and three move off the imaginary axis The analysis suggests that the energy 
equation with enthalpy damping be written in the form 

(pE)i + (puH)j + (pvH)y + ap(E + kc 2 /~f)(H - Hoo) = 0 
The eigenvalues that move off the imaginary axis move into the left half-plane provided that k > 1 and u 2 + v 2 < c 2 
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It is appropriate here to comment on the missing factor in the analysis, the dissipation term Going back to (1.7), 
take a dissipation term of the form 

D{w) = (d* + d*)(w u w 3 ,Wi,pH) T 

The point here is that the fourth component is not 1^4, it is instead pH. This is done because the system of equations 
with artificial dissipation still admits steady solutions with H = Hqq. To see this, note that when the continuity 
equation with dissipation is multiplied by Hoc and subtracted from the energy equation with dissipation we obtain 

( P (E - J7oo)) f + (Htf “ tfoo))* + (PV{H - Hoofty 

+ ap{h A - tfoo)(tf - ffoo) + t[d\ + dl)p{H - Hoc) = 0 

which evidently admits steady solutions with H = H^. We would want to look at the polynomial 

p(z) .= d et(zl + i(uiA + u 2 B) + aC + + t^)^) 

where D = dD/dxv. If we do this and then form N~ l DN y we find 


N~ l DN 


/ 1 0 0 
0 10 
0 0 1 
\ (7 — l) 2 (u 2 + v 2 )/2 0 0 



The fact that the fourth row here is not the fourth row of the identity matrix has the consequence that p(z) with 
e ^ 0 is vastly more complicated than p(z) with e = 0 This author has not been able to make any headway in 
analyzing the case f / 0 because of the formidable algebraic difficulties In fact, even for the case of no enthalpy 
damping (a = 0) this author is not aware of any proof of stability of the linearized system with dissipation, because 
the dissipation term is applied to a nonlinear function of w 

Some experiments were done with the code FL052R in changing the form of the enthalpy damping term in the 
energy equation The form hi = E -f kc 2 / 7 was used, and various values of k were tried Only one line of the code 
was changed The first case was flow about a NACA0012 airfoil at a freestream Mach number 0 8 and at an angle of 
attack of 1 25° The code was run in multigrid mode with 40 cycles on a 48 X 8 mesh, then 40 cycles on a 96 X 16 
mesh, and finally 100 cycles on a 192 x 32 mesh. A five-stage scheme with coefficients (l/4, 1/6, 3/8, 1/2, 1) was used, 
with two dissipation evaluations per time step For the unmodified code, the convergence rate for 100 cycles on the 
192 x 32 grid was 0 9557 A table of convergence rate vs k for the modified code is shown here 

k 1 00 1 25 1 50 1 75 2 00 2 50 3 00 4 00 

convergence rate . 0 9864 0 9565 0 9538 0 9553 0 9580 0 9630 0 9672 0 9788 

Evidently the optimum k is approximately 1 5 For a subcntical case with a freestream Mach number of 0 63 and 
an angle of attack of 2°, with all other parameters the same, the unmodified code had convergence rate 0 9288 The 
modified code produced the following results' 


k 1.00 1 50 2 00 2 50 3 00 3 50 4 00 

convergence rate . . 0 9451 0.9303 0 9282 0 9276 0 9273 0.9275 0 9302 


Here the optimal k is apparently between 2 5 and 3 5 In both cases it is clearly difficult to improve on the original 
code This author believes the excellent convergence behavior of the unmodified code may be explained as follows. 
The energy equation enthalpy damping term in the unmodified code comes from discretizing 


d t (p(E - /Too)) + d x {pn(H - Hoo)) + d y {pv{H - H „)) + ap{H - H^) = 0 

1 e , this equation is used in place of the fourth equation of (2 9). In terms of the partial differential equation system, 
an equivalent equation is obtained by multiplying the first equation of (2 9) by and adding to this, obtaining 

dt(pE) -f d x (puH) + d y {pvH) + ap(H - i7oo)(i?oo + 1) = 0 


9 



In the notation of (3 2), h 4 = Hoo + 1. With the definition of /j 4 as h 4 = E + kc 2 / 7 we would obtain h 4 = H 00 + 1 
if we had k = 7 (jEoo — i? + c ^/7 + 1 ) / c 2 . Now, the nondimensionalization in the FL052 codes is such that = 7 , 
hence at freestream conditions we would have k = 2 . Thus h 4 = E + 2^/7 is roughly the same as the unmodified 
code; also the previously noted experiments showed k = 1 5 in one case and k G (2 5, 3 5) in the other case were 
near-optimal. Thus k = 2 (nearly the original code) would be expected to give excellent results It is also clear 
that the unmodified code is equivalent to a code with locally varying k ; this author has not experimented with other 
forms of locally varying enthalpy damping terms. 


4. Finite Difference Formulation with Enthalpy Damping 

The aim of this section is to study the actual finite difference implementation of enthalpy damping. The aim is to 
understand why enthalpy damping can be effective in practice when on the level of the partial differential equation 
it does not produce a system with solutions which tend to zero. 

The general formulation of enthalpy damping to be studied is as follows. First spatial differencing is applied, 
sending tv n to uJ" := 7(w n ). Then enthalpy damping is applied, via 

u / n+1 - w” + aAfJ/(u; n+1 , w*) = 0 ( 4 . 1 ) 

For example, the first equation could be taken to be 

P" +1 -P?,+ «A <P" +1 (^ - tfoo) = 0 

Notice that in this formulation (the FL052 formulation) the enthalpy damping term is in a semi-implicit form which 
is readily solved for p* +1 . If the enthalpy damping term were fully implicit, it would be nonlinear in w n+l and 
difficult to solve 

We plan to study the stability of (4 1 ) about a solution w* with tv * = 7(u;*) and #(«;*,«;*) = 0 . Perturbing to 
w* + eu;, we find the perturbation equation is 

ti > ft+1 - Atv n + aAt(Biiv n+l *f BiAtv") =0 (4 2 ) 


where 

A = DT(tv*) (the Jacobian matrix) 

B\ = Di)l(w* ^w*) (Jacobian with respect to first argument) 

= Djl/fu;*, u/*) (Jacobian with respect to second argument) 

The equation (4 2) is equivalent to 

(J + aAfBi)u > n+1 = (J - a&t 82 )Aw n (4.3) 


Let us take 

#K +1 ,^) = {H» - tfoo)*'K + 1 ,t^) 

Then B\ = 0 and Bj = J/'(tu*, u;*) • V~( H n — # 00 )* Furthermore, if the basic iteration scheme is an explicit 
multistage method with central spatial differencing, then T(tv n ) = P(AtL)w n , where P is a polynomial, L = 
A8 x j2Ax + B6 y /2Ay, and 8 X and S v are central space difference operators Thus the matrix problem (4 3) becomes 

tr;"* 1 =(I-aAtBi)P(AtL)w n 

Fourier transforming in space produces the problem 

u ) n+1 = (7 - a&tB 2 )P{AtL)w n 

where A IL = i(uiA + LJ 2 B) and = At sin 6/ Aar, cj 2 = At sin0/Ay, and 0 , <j> are the Fourier dual variables. 
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Figure 4 1. Contours of |P(^)| in the complex plane for classical 
fourth-order Runge-Kutta. Contour levels: 0 1 , 0.2, 0 3,..., 10 . 

A reasonable form (consistent with the work of section 3) for the implementation of the enthalpy damping step is 

p n+1 -fi + aAtp n + 1 (HZ-H co )= 0 

(pu) n+1 - (pu) n + aAt(pu) n+l (H " - tfoo) = 0 

/ - Sw/ l * *) 

(pt>) n+1 - (pv) n + aA<(pt/) n+1 (.ff n - /Too) = 0 
(pE) n+1 - (rf) n + aAf[(p£) n+I + *p*](E* -Hoo)=0 
The parameter k here is the same as in the previous section Then our perturbation equation is 

u? n+1 = (I -aAtBi)P(AtL) 

For the decay of perturbations about tr*, it is necessary and sufficient that all eigenvalues A of (7 — aAtB 2 )P[AtL) 
satisfy |A| < 1 Furthermore, the asymptotic rate of decay is given by the modulus of the largest eigenvalue. For 
a = 0 (no enthalpy damping) we recover the usual CFL-like condition max^ |P(Aj)| < 1, where A^ are the eigen\alues 
of AtL : Ai = Aj = ivU, Aj = iu(U + c), A4 = iu[U — c), where uU := u?i u + u 2 v and u/ 2 := u\ + u\. 

For example, the explicit Euler method has P(AtL) = I — AtL (this is unstable); the classical fourth-order 
Runge-Kutta method has P{x) = 1 - x + x 2 j2 - z 3 /6 + x*/24, for which it is known that the stability condition 
max, |^(Aj)| < 1 becomes max, |Aj| < 2y/2 . Contours of \P(z)\ in the complex z-plane are shown in Figure 4 1 for 
the fourth-order Runge-Kutta method. 
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The question before us now is, does the modulus of the maximum eigenvalue decrease when a At is increased from 
0? Also, what are good values for a At and kl For U as in ( 4 . 4 ) we have 




( P \ 
pu 

pv 

\pE + kpJ 


. i(_ 7 £ + (7 - l)(u 2 + v 2 ), (1 - 7)u, (1 - 7)1;, 7) 


Again we have a matrix of rank one, B2, and a diagonalizable matrix P(AtL) with a double eigenvalue at P(iuU), 
so the theorem of the previous section implies that (7 — aAtB2)P(AtL) has as one of its eigenvalues P(twU). Thus 
if this eigenvalue is the eigenvalue of maximum modulus, then enthalpy damping will not improve the convergence 
rate of the iteration. On the other hand, if one of the other eigenvalues is obstructing the convergence, then enthalpy 
damping may help. 

To help in studying this question, it is useful to diagonalize the matrix L. The matrix T which diagonalizes 
ui A 4- u; 2 B is known analytically (see Warming and Beam (Ref 14 )). It is 



f 1 0 

p/(y/ 2 c) 

p/s/ 2 c ^ 

T — 

0 W2/W 

UiJ{s/ 2 u) 

— Wi/(s/ 2 u) 

1 — 

0 -Ui/u 

w 2/(v^ w ) 

— U 2 /(s/ 2 u) 


^0 0 

pcf s /2 

Pc/sj 2 j 

T~ l (I - aAtBi)P(AtL)T = (I - aAtT~ l B^T)T' 

~ l P(&tL)T 


= (7 — aAtT~ x B^T) 

• dia g{/>(iW), />(*W), P{iu[U 4- c]), P(iu[U - c})}. 

The characteristic polynomial of this matrix vanishes at P(tuU), and upon division by z - P(iuU) the characteristic 
polynomial becomes (writing a = P(iuU), <7+ = P(\u\U + c]), a _ = P{xu\U - c]), and M = 77 /c), 


(z - a)(z - <T+){z - (T-) 

/y /\ 

4 [z 2 {2(k - l)<7 -f (1 4- M)(ik - k + 1)<t+ 4- (l - M)(^k - k 4- l)<r_}- 

27 

- ^({7^(1 + M) 4 - (fc - 1)(1 - A/)}<r<7+ 4- {7^(7 “ M) + {k - 1)(1 4- \I)}(T<r _ 

4* 2(7 k - k + ) 

4" 27fc<7<7+<T_] 


(4 5 ) 


Note that a has units of time _1 veIocity“ 2 (cf. (4 4 )), so a Ate 2 is dimensionless. 

We will study the roots of the polynomial in (4 5 ) numerically. Fix a “Mach number” A/ = U/c. Define v (“CFL 
number”) as max(|wf/|, \u(U 4* c)|, \u(U - c)|). Then 

uU= /^sgn(A/)/max(|l + l/M|,|l-l/M|), ifA/#0 f46) 

1 0 if M = 0 ' ’ 

Now for a given M we can graph the spectral radius vs u by using the definition of uU from (4 6) and then numerically 
finding the roots of the polynomial (4 5 ). We can further search for an optimal value of the dimensionless parameter 
a Ate 2 by using as optimality criterion the integral between the curves of spectral radius vs. v for a Ate 2 = 0 and 
a Ate 2 ^ 0, and searching for a Ate 3 which maximizes this integral (Numerically, the composite trapezoidal rule was 
used for the integration.) Two typical graphs of spectral radius vs. v without and with enthalpy damping are shown 
in Figure 4 2 Note that enthalpy damping has little or no effect at small CFL numbers, it is only for CFL numbers 
near the stability limit that enthalpy damping is significantly helpful For supersonic flow (M > 1), the spectral 
radius with enthalpy damping was observed to be greater than 1 . In all subsonic cases investigated, the spectral 
radius with enthalpy damping was exactly the same as or slightly less than the spectral radius without enthalpy 
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CFL NUMBER 

Figure 4 2(a). Spectral radius vs. CFL number for numerical scheme without and with enthalpy- 
damping. Five-stage method with coefficients (1/4, 1/6, 3/8, 1/2, l),Af = 0 8, k = 2, aAtc 2 = 0 25 
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Figure 4 2(b). Same as (a) except aA tc 2 =0 5. 

damping, at least in the region of stability The five-stage method with coefficients (1/4, 1/6, 3/8, 1/2, 1) was studied, 
with 0 < v < 4 (the stability limit) With k = 2 and two different Mach numbers, A/ = 0 8 and A/ = 0 4, the 
following results were obtained for the area between the curves as a function of a Ate" 2 . 

aAte 2 0 30 0 40 0.50 0 60 0.70 0 80 0 90 100 

area (A/ = 0 8) . . . 0 0840 0.0967 0.101 0103 0104 0105 0.105 0 0685 

area (M = 0.4) . . . 0.0039 0 0039 0 0039 0 0039 0 0039 0 0039 0 0039 0.00277 
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Evidently there is not a sharp optimum value of a Ate 2 , for in the first case any a Ate 2 between 0 5 and 0 9 gives 
approximately the same result, while in the second case any a Ate 2 E (0.3, 0.9) can be chosen. If the speed of sound 
is normalized to be 7 at infinity (as in the FL052 codes), then our good range for a At becomes a At E (0 36, 0 64) 
in the first case, and a At E (0.21,0.64) in the second case. 

5. Conclusions 

The Euler equations with enthalpy damping can be studied as a system of partial differential equations. Linear 
stability analysis of the system around a state with constant enthalpy shows that not all solutions decay to zero. The 
analysis suggests that the enthalpy damping term in the energy equation be written as ap[E + kc 2 /i)(H - Hoo). 
For this formulation, the linear stability analysis shows stability (neutral stability, i e., no eigenvalues in the open 
right half plane) if k > 1 and the flow is subsonic. 

When the equations are discretized via an explicit multistage method, the spectral radius can be decreased by 
enthalpy damping, even though it is again the case that one of the eigenvalues is unchanged. A numerical study with 
k = 2 revealed that a Ate 2 « 0 5 should give good results The particular value of a At this leads to depends on the 
nondimensionahzation of a given computer code 

This analysis has not taken the dissipation term into account, because the dissipation operator m the energy 
equation is applied to pH , not pE (so that the equations with dissipation will still have constant enthalpy at 
steady state), which makes the analysis much more complicated All the analysis shows that enthalpy damping 
should be destabilizing for supersonic flow, but in practice enthalpy damping does not harm the convergence for the 
FL052 codes in the transonic range. It is probably the case that the dissipation term is effective in stabilizing the 
computation, even in locally supersonic regions of flow 
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